Generate data

# generate some triangle-shaped data
set.seed(234532)
np$random$seed(134223L)
# generate archetypes
arc = np$column_stack(list(np$random$normal(20, 5, 4L),
                           np$random$normal(60, 5, 4L),
                           np$random$normal(40, 5, 4L)))
arc = apply(arc, 1, sample)
# generate weights
sim_weights = np$random$dirichlet(c(0.1, 0.1, 0.1), 10000L)
data = np$matmul(sim_weights, arc)
data = np$random$poisson(data)
rownames(data) = paste0("C", 1:nrow(data))
data = t(data)
dim(data)
## [1]     4 10000

Standard Poisson PAA

# Poisson archetypal analysis
arc_data = fit_pch(data, 3L, method = "poisson_aa",
                   method_options = list(model_fun = ParetoTI::paa_poisson,
                                         c_alpha_prior = 0.0001,
                                         weight_alpha_prior = 0.1,
                                         #greta::adam(learning_rate = 0.1)),  greta::l_bfgs_b(),  greta::momentum(learning_rate = 0.001, momentum = 0.9, use_nesterov = TRUE)
                                         optimiser = greta::adam(learning_rate = 0.1),
                                         # initialise archetypes at louvain cluster centers
                                         initial_values = "louvain_centers"
                                         #covar = covar
                   ),
                   conv_crit = 1e-04,
                   maxiter = 5000, volume_ratio = "none")
## Registered S3 method overwritten by 'R.oo':
##   method        from       
##   throw.default R.methodsS3
plot_arc(arc_data, data, data_size = 0.1, geom = geom_bin2d)

Poisson PAA with no archetype constraint

# Poisson archetypal analysis with no archetype position constraint
arc_data_free = fit_pch(data, 3L, method = "poisson_aa",
                        method_options = list(model_fun = ParetoTI::paa_poisson_free,
                                              weight_alpha_prior = 0.01,
                                              #greta::adam(learning_rate = 0.1)),  greta::l_bfgs_b(),  greta::momentum(learning_rate = 0.001, momentum = 0.9, use_nesterov = TRUE)
                                              optimiser = greta::adam(learning_rate = 0.1),
                                              # initialise archetypes at louvain cluster centers
                                              initial_values = "louvain_centers"
                                              #covar = covar
                        ),
                        conv_crit = 1e-04,
                        maxiter = 5000, volume_ratio = "none")

plot_arc(arc_data_free, data, data_size = 0.1, geom = geom_bin2d)

PCHA

# PCHA for normal data
arc_data_pcha = fit_pch(data, 3L, method = "pcha",
                        conv_crit = 1e-06,
                        maxiter = 5000, volume_ratio = "none")
plot_arc(arc_data_pcha, data, data_size = 0.1, geom = geom_bin2d)

MCMC for poisson paa

m = paa_poisson_free(t(data), n_arc = 3, weight_alpha_prior = 0.8, scale_data_sd = 0.2)
#m = paa_poisson(t(data), n_arc = 3, weight_alpha_prior = 0.8, c_alpha_prior = 0.01)
draws <- mcmc(m$model, n_samples = 1000, chains = 4)

#draws_archetypes = draws[,grepl("archetypes", colnames(draws[[1]]))]
draws_archetypes = calculate(m$archetypes, draws)
bayesplot::mcmc_trace(draws_archetypes)

# extract samples in the standardised form
library(foreach)
reshape_draws = function(parameter = "archetypes", draws, m){
  draws_param = calculate(m[[parameter]], draws)
  # combine chains
  draws_param = foreach(i = seq_along(draws_param), .combine = rbind) %do% draws_param[[i]]
  # reshape into dim(chains * samples, param_dim_1, param_dim_2)
  draws_param = array(draws_param, dim = c(nrow(draws_param),
                                           dim(m[[parameter]])[1],
                                           dim(m[[parameter]])[2]))
  # convert array to list
  foreach(i = seq_len(dim(draws_param)[1]), .combine = c) %do% list(t(draws_param[i,,]))
}
draws_archetypes = reshape_draws("archetypes", draws, m)
draws_archetypes = lapply(draws_archetypes, exp)
draws_weights = reshape_draws("weights", draws, m)
#draws_c = reshape_draws("c", draws, m)

# align to reference - first polytope from bootstraping
ref_XC = draws_archetypes[[1]]
for (i in seq_along(draws_archetypes)) {
  ind = align_arc(ref_XC, draws_archetypes[[i]])$ind
  draws_archetypes[[i]] = draws_archetypes[[i]][, ind]
  draws_weights[[i]] = draws_weights[[i]][ind, ]
  #draws_c[[i]] = draws_c[[i]][ind, ]
}

# PCHA with bootstrapping
arc_data_pcha = fit_pch_bootstrap(data, n=3,
                                  noc = 3L, method = "pcha",
                                  conv_crit = 1e-06,
                                  maxiter = 5000, volume_ratio = "none")
plot_arc(arc_data_pcha, data, data_size = 0.1)

# replace pcha results with MCMC draws
arc_data_pcha$pch_fits$XC = draws_archetypes
arc_data_pcha$pch_fits$S = draws_weights
#arc_data_pcha$pch_fits$C = draws_c

plot_arc(arc_data_pcha, data, data_size = 0.1, arch_size = 0.1, line_size = 0)

Normal PAA with no archetype constraint

####### ---- Normal model ---- #######
m = paa_normal_free(t(data),
                    n_arc = 3,              # number of achetypes
                    weight_alpha_prior = 0.8,
                    scale_data_sd = 0.2,
                    covar = NULL)
m$arc = arc
plot(m$model)
####### ---- finding parameters by optimisation ---- #######
set.seed(100)
# optimise instead of MCMC
time = Sys.time()
evalq({
  opt_res = opt(model,
                optimiser = adam(learning_rate = 0.1), #adadelta(learning_rate = 0.001) adam(learning_rate = 0.001), l_bfgs_b()   # optimisation method used to find prior-adjusted maximum likelihood estimate: adam, l_bfgs_b
                initial_values = initials(archetypes_mean = arc),
                max_iterations = 10000,
                tolerance = 1e-4, adjust = TRUE,
                hessian = FALSE)
}, envir = m)
Sys.time() - time
## Time difference of 1.85192 mins
opt_res = m$opt_res
opt_res$iterations
## [1] 10000
arc_norm_free = copy(arc_data_free)
arc_norm_free$XC = t(opt_res$par$archetypes_mean)

# Generate data from the model parameters
evalq({
  mu_sample = calculate(mu, values = opt_res$par)
  sd_sample = calculate(sd, values = opt_res$par)
}, envir = m)
data_sample = array(rnorm(length(m$mu_sample), m$mu_sample, m$sd_sample), dim = dim(m$mu_sample))

cowplot::plot_grid(
  plot_arc(arc_norm_free, data, data_size = 0.1,
           geom = ggplot2::geom_bin2d) +
    xlim(-10, 140) + ylim(-10, 100) +
    ggtitle("Data"),
  plot_arc(arc_norm_free, t(data_sample), data_size = 0.1,
           geom = ggplot2::geom_bin2d) +
    xlim(-10, 140) + ylim(-10, 100) +
    ggtitle("Sample from the model"),
  ncol = 2
)
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_text).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_text).

Negative binomial PAA with no archetype constraint

####### ---- Negative binomial model ---- #######
m = paa_nb_free(t(data),
                n_arc = 3,              # number of achetypes
                weight_alpha_prior = 0.8,
                scale_data_sd = 0.2,
                covar = NULL)
m$arc = arc
plot(m$model)
####### ---- finding parameters by optimisation ---- #######
set.seed(100)
# optimise instead of MCMC
time = Sys.time()
evalq({
  opt_res = opt(model,
                optimiser = adam(learning_rate = 0.1), #adadelta(learning_rate = 0.001) adam(learning_rate = 0.001), l_bfgs_b()   # optimisation method used to find prior-adjusted maximum likelihood estimate: adam, l_bfgs_b
                #initial_values = initials(archetypes_mean = arc),
                max_iterations = 10000,
                tolerance = 1e-4, adjust = TRUE,
                hessian = FALSE)
}, envir = m)
Sys.time() - time
## Time difference of 2.038677 mins
opt_res = m$opt_res
opt_res$iterations
## [1] 7206
arc_norm_free = copy(arc_data_free)
arc_norm_free$XC = t(opt_res$par$archetypes_mean)

evalq({
  size_sample = calculate(size, values = opt_res$par)
  prob_sample = calculate(prob, values = opt_res$par)
}, envir = m)
data_sample = array(rnbinom(length(m$size_sample), size = m$size_sample, prob = m$prob_sample), dim = dim(m$size_sample))

cowplot::plot_grid(
  plot_arc(arc_norm_free, data, data_size = 0.1,
           geom = ggplot2::geom_bin2d) +
    xlim(-10, 160) + ylim(-10, 100) +
    ggtitle("Data"),
  plot_arc(arc_norm_free, t(data_sample), data_size = 0.1,
           geom = ggplot2::geom_bin2d) +
    xlim(-10, 160) + ylim(-10, 100) +
    ggtitle("Sample from the model"),
  ncol = 2
)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 11 rows containing non-finite values (stat_bin2d).
## Warning: Removed 1 rows containing missing values (geom_tile).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).

Poisson PAA with no archetype constraint (posterior predictive check)

####### ---- Poisson model ---- #######
m = paa_poisson_free(t(data),
                     n_arc = 3,              # number of achetypes
                     weight_alpha_prior = 0.8,
                     scale_data_sd = 0.2,
                     covar = NULL)
m$arc = arc
plot(m$model)
####### ---- finding parameters by optimisation ---- #######
set.seed(100)
# optimise instead of MCMC
time = Sys.time()
evalq({
  opt_res = opt(model,
                optimiser = adam(learning_rate = 0.1), #adadelta(learning_rate = 0.001) adam(learning_rate = 0.001), l_bfgs_b()   # optimisation method used to find prior-adjusted maximum likelihood estimate: adam, l_bfgs_b
                #initial_values = initials(archetypes = arc),
                max_iterations = 5000,
                tolerance = 1e-4, adjust = TRUE,
                hessian = FALSE)
}, envir = m)
Sys.time() - time
## Time difference of 54.37565 secs
opt_res = m$opt_res
opt_res$iterations
## [1] 5000
arc_norm_free = copy(arc_data_free)
arc_norm_free$XC = t(opt_res$par$archetypes)

evalq({
  mu_sample = calculate(mu, values = opt_res$par)
}, envir = m)
data_sample = array(rpois(length(m$mu_sample), lambda = m$mu_sample), dim = dim(m$mu_sample))

cowplot::plot_grid(
  plot_arc(arc_norm_free, data, data_size = 0.1,
           geom = ggplot2::geom_bin2d) +
    xlim(-10, 140) + ylim(-10, 100) +
    ggtitle("Data"),
  plot_arc(arc_norm_free, t(data_sample), data_size = 0.1,
           geom = ggplot2::geom_bin2d) +
    xlim(-10, 140) + ylim(-10, 100) +
    ggtitle("Sample from the model"),
  ncol = 2
)